

\chapter{偏微分方程}

我们依照\cite{Hjorth-Jensen}，简单介绍扩散方程、拉普拉斯方程、波动方程的有限差分解法。

\section{扩散方程}
在数学物理方法的课程中，我们学过热量扩散方程,
\begin{equation}
\nabla^2 T  = D \frac{ \partial T }{ \partial t },
\end{equation}
如果稍微重新定义$T$，即可把微分方程写成如下形式,
\begin{equation}
\nabla^2 u = \frac{ \partial u }{ \partial t }，
\end{equation}
一维情况下，扩散方程为，
\begin{equation}
\frac{ \partial^2 u }{ \partial x^2 } = D \frac{ \partial u }{ \partial t },
\end{equation}
$u(x,t)$在$x$轴上的取值范围是$[x_a,x_b]$，在$t$轴上的取值范围是$[t_a,t_b]$，在初始时刻,
\begin{equation}
u(x,t=t_a) = g(x),
\end{equation}
在边界上，
\begin{eqnarray}
u(x_a,t) &=& a(t), \\
u(x_b,t) &=& b(t).
\end{eqnarray}
这样就限定了唯一的解。

\paragraph{显式递推公式}
利用有限差分，可以近似处理偏导数，
\begin{eqnarray}
\frac{ \partial^2 u }{ \partial x^2 } &=& \frac{1}{(\Delta x)^2}( u(x+ \Delta x,t) + u(x- \Delta x,t) - 2u(x,t) ) + O((\Delta x)^2),
\\
\frac{ \partial u }{ \partial t } &=& \frac{1}{\Delta t}( u(x,t+ \Delta t) - u(x,t) ) + O(\Delta t).
\end{eqnarray}
带入一维扩散方程，得到
\begin{equation}
\frac{1}{(\Delta x)^2}( u(x+ \Delta x,t) + u(x- \Delta x,t) - 2u(x,t) )
\approx \frac{1}{\Delta t}( u(x,t+ \Delta t) - u(x,t) ),
\end{equation}
误差为$O((\Delta x)^2) + O(\Delta t)$。
记$\alpha = \frac{ \Delta t }{ (\Delta x)^2 }$，则方程变为
\begin{equation}
u(x,t+ \Delta t) \approx \alpha u(x+ \Delta x,t) + (1- 2\alpha)u(x,t)+ \alpha u(x- \Delta x,t),
\end{equation}
这意味着，如果我们知道了$t$时刻的情况，就可以往前推断$\Delta t$之后的情况。
若$u(x,t)$在$x$轴上的取值范围是$[x_a,x_b]$，在$t$轴上的取值范围是$[t_a,t_b]$，我们把这个定义域用四方网格划分，格点在$x$方向上间距为$\Delta x$，在$t$方向上的间距为$\Delta t$，记
\begin{eqnarray}
x_i &=& x_a + i*\Delta x, \\
t_j &=& t_a + j*\Delta t,
\end{eqnarray}
格点上的$u$值标记为
\begin{equation}
u_{i,j} = u(x_i, t_j) \approx u(x_a + i*\Delta x, t_a + j * \Delta t ),
\end{equation}
那么，数值递推式为
\begin{equation}
u_{i,j+1} = \alpha u_{i-1,j} + (1- 2\alpha)u_{i,j} + \alpha u_{i+1,j},
\end{equation}
根据前人的推断\cite{Golub2013Matrix}，这个公式在如下条件下是稳定的，
\begin{equation}
\Delta t / (\Delta x)^2 \leq 1/2,
\end{equation}
所以如果$\Delta t$太大，结果可能是错误的。

\paragraph{隐式递推公式}
在上面的推导中，我们计算$u$对时间的偏导数时，用了公式
\begin{equation}
\frac{ \partial u }{ \partial t } = \frac{1}{\Delta t}( u(x,t+ \Delta t) - u(x,t) ) + O(\Delta t).
\end{equation}
我们也可以用这个公式
\begin{equation}
\frac{ \partial u }{ \partial t } = \frac{1}{\Delta t}( u(x,t) - u(x,t- \Delta t) ) + O(\Delta t).
\end{equation}
那么，经过类似的推导，并整理后，得到数值递推式为
\begin{equation}
u_{i,j-1} = -\alpha u_{i-1,j} + (1+ 2\alpha)u_{i,j} - \alpha u_{i+1,j},
\end{equation}
如果$u_{*,j-1}$都已经已知，我们可以得到关于$u_{*,j}$的$N_x$个线性方程，$N_x$为$x$轴上$(x_a,x_b)$之间的格点数。
\begin{equation}
\begin{bmatrix}
1+2\alpha & -\alpha \\
-\alpha   & 1+2\alpha & -\alpha \\
          & -\alpha   & 1+2\alpha & -\alpha \\
          &			  &	\ddots	  & \ddots  & \ddots \\
          &			  &			  & -\alpha & 1+2\alpha
\end{bmatrix}
\begin{bmatrix}
u_{1,j} \\
u_{2,j} \\
\vdots \\
u_{N_x,j}
\end{bmatrix}
=
\begin{bmatrix}
u_{1,j-1} + \alpha a(t_j)\\
u_{2,j-1} \\
\vdots\\
u_{N_x,j-1} + \alpha b(t_j)
\end{bmatrix}
\end{equation}
于是我们用三对角线性方程组的数值解法，求解上面的方程组，就得到了$u_{*,j}$，因为上式没有写出$u_{*,j}$的显式表达式，所以这种方法叫做隐式递推。
前文中提到过，解三对角线性方程组只需要$O(N_x)$次操作，所以是比较有效的。
根据前人的证明，这个算法是稳定的。

下面是隐式递推的代码，征用了前面章节编写的三对角方程组求解函数TridiagLinear。
\begin{lstlisting}
/*
 * diffusion1d(...)             solves a 1-dimensional diffusion equation
 *              \partial^2 u / \partial x^2 = \partial u / \partial t
 *              and prints u(*,t_j) into "diffusion1d.txt"
 * double xa
 * double xb                    x zone [xa, xb]
 * double (*g_initial)(...)     u(x, t=0)
 * double ta
 * double tb                    t zone [ta, tb]
 * double (*boundary_a)(...)    u(xa, t)
 * double (*boundary_b)(...)    u(xb, t)
 * int Nx                       grid point number in (xa,xb)
 * int Nt                       grid point number in (ta,tb)
 */
void diffusion1d( double xa, double xb, double (*g_initial)(double x), double ta, double tb, double (*boundary_a)(double t), double (*boundary_b)(double t), int Nx, int Nt){

        int i, j;
        double dx, dt, alpha;
        double *u = new double [Nx];
        double *unext = new double [Nt];
        double *a = new double [Nx];// tridiagonal matrix
        double *b = new double [Nx];
        double *c = new double [Nx];

        dx = (xb-xa)/(Nx+1);// [xa, xb] is divided into Nx+1 small zones
        dt = (tb-ta)/Nt;// [ta, tb] is divided into Nt steps
        alpha = dt/dx/dx;
        for(i=0;i<Nx;i++){// initialize tridiagonal matrix
                a[i] = -alpha;
                b[i] = 1 + 2*alpha;
                c[i] = -alpha;
        }
        for(i=0;i<Nx;i++){//initial u(x, t=0)
                u[i] = g_initial(xa + (i+1)*dx);
        }
        ofstream fp("diffusion1d.txt");
        for(j=1;j<=Nt;j++){
                fp<< boundary_a( ta+ (j-1)*dt)<<"\t";//output u(*,j-1) into file
                for(i=0;i<Nx;i++)fp<<u[i]<<"\t";
                fp<< boundary_b( ta+ (j-1)*dt)<<"\t";
                fp<<endl;
                u[0] += alpha * boundary_a( ta + j*dt);
                u[Nx-1] += alpha * boundary_b( ta + j*dt);
                TridiagLinear( Nx, a, b, c, unext, u);
                for(i=0;i<Nx;i++) u[i] = unext[i];
        }
        fp.close();
        delete [] c;
        delete [] b;
        delete [] a;
        delete [] u;
        delete [] unext;
}
\end{lstlisting}

自定义初始条件，边界条件以后，就可以调用上面的代码计算1维扩散问题。
下面的代码做一个热量传导问题，
\begin{lstlisting}
头文件。。。
函数diffusion1d(...)

double Tleft(double t){// temperature of left boundary
        return 100;
}

double Tright(double t){// temperature of right boundary
        return 0;
}

double Tinitial(double x){
        return 0;
}

int main(){
        diffusion1d(0, 1, Tinitial, 0, 0.5, Tleft, Tright, 1000, 1000);
        return 0;
}
\end{lstlisting}
冬天非常寒冷，一根铁棒上各处温度都是0摄氏度，右端浸泡在0摄氏度的水里，左端与100摄氏度的沸水接触，随着时间的过去，热量会源源不断地从铁棒左端流向右端，铁棒上的温度分布会趋于均匀。

运行上面的代码以后，生成diffusion1d.txt文件，然后用gnuplot，可以画出热量图，如图\ref{fig:diffusion1d}
\begin{figure}[h]
\centering
\includegraphics[width=0.8\textwidth]{diffusion1d}
\caption{一维热扩散问题：铁棒上的温度演化图。铁棒左端恒为100摄氏度，右端恒为0摄氏度，初始时铁棒上温度为0。}
\label{fig:diffusion1d}
\end{figure}

\section{二维泊松方程}
二维泊松方程为
\begin{equation}
\nabla^2 u = \rho,
\end{equation}
静电势的方程就是这个形式。
这里我们求解二维情形，三维求解的道理也一样，只是计算量比较大，算一次要等很久。

二维情况下，泊松方程变为
\begin{equation}
\frac{\partial^2 }{\partial x} u + \frac{\partial^2}{\partial y} u = \rho.
\end{equation}
用有限差分公式代替二阶偏导数，得到
\begin{equation}
\frac{ u_{i+1,j} + u_{i-1,j} - 2 u_{i,j} }{(\Delta x)^2} + O((\Delta x)^2)
+ \frac{ u_{i,j+1} + u_{i,j-1} - 2 u_{i,j} }{(\Delta y)^2} + O((\Delta y)^2) = \rho_{i,j}
\end{equation}
其中，$u_{i,j} = u(x_a + i\Delta x, y_a + j\Delta y)$.
那么，如果取$\Delta x = \Delta y = h$，在2阶精度下，数值递推式为
\begin{equation}
u_{i-1,j} + u_{i+1,j} + u_{i,j+1}, u_{i,j-1} - 4 u_{i,j} = h^2 \rho_{i,j}.
\end{equation}
写成矩阵会得到一个关于$u_{i,j}$的线性方程组，系数矩阵是一个稀疏矩阵，求得线性方程组的解，就得到任意格点的$u_{i,j}$，就认为得到了问题的数值解。

如果用高斯消元法，或者LU分解，会得到大量的非零矩阵元，使整个问题变得更加复杂。
所以这种问题，更适合用迭代法。
如果用Jacobi方法，即
\begin{eqnarray}
A \vec{x} &=& \vec{b}, \\
A &=& U + D + L, \\
\vec{x}_{j+1} &=& D^{-1}( \vec{b} - (U+L) \vec{x}_j ).
\end{eqnarray}

给每个格点上的$u_{i,j}$任意赋一个初始值，然后用Jacobi方法进行迭代，实际上就相当于
\begin{equation}
u^{(k+1)}_{i,j} = \frac{1}{4}(u^{(k)}_{i-1,j} + u^{(k)}_{i+1,j} + u^{(k)}_{i,j+1}, u^{(k)}_{i,j-1} - h^2 \rho_{i,j}).
\end{equation}
如果在计算等号左边的$u^{(k+1)}_{i,j}$时，右边的式子中有些项已经做了第$k+1$轮迭代，比如，已经有$u^{(k+1)}_{i-1,j}$，那么可以不用$u^{(k)}_{i-1,j}$，而用$u^{(k+1)}_{i-1,j}$，这样做就是高斯-席德尔方法。

我们可以用高斯席德尔方法，写一个通用的二维四方网格扩散方程求解函数，然后做一个简单的例子：匀强电场。
\begin{lstlisting}
void poisson2d(int n, double **u, double xa, double ya, double h, double precision,
                double (* boundary_x1)(double x), double (* boundary_x2)(double x),
                double (* boundary_y1)(double y), double (* boundary_y2)(double y),
                double (* rho)(double x, double y) ){

        int i,j;
        double temp, change = precision +1;

        for(i=0;i<n;i++) u[i][0] = boundary_x1( xa + i*h );
        for(i=0;i<n;i++) u[i][n-1] = boundary_x2( xa + i*h );
        for(j=0;j<n;j++) u[0][j] = boundary_y1( ya + i*h );
        for(j=0;j<n;j++) u[n-1][j] = boundary_y2( ya + i*h );

        for(i=1;i<n-1;i++)
                for(j=1;j<n-1;j++)
                        u[i][j]=0;

        while( change > precision ){

                change = 0;
                for(i=1;i<n-1;i++){
                        for(j=1;j<n-1;j++){
                                temp = u[i][j];
                                u[i][j] = 0.25 * (u[i-1][j] + u[i+1][j] + u[i][j-1] + u[i][j+1] - h * h * rho(xa+i*h, ya+i*h) );
                                change += fabs( u[i][j] - temp );
                        }
                }
        }
        ofstream fp("poisson2d.txt");
        for(j=0;j<n;j++){
                for(i=0;i<n;i++){
                        fp<<u[i][j]<<"\t";
                }
                fp<<endl;
        }
        fp.close();
}
\end{lstlisting}
要做匀强电场，需要设置恰当的边界条件，
\begin{lstlisting}
double x1(double x){
        return x;
}
double x2(double x){
        return x;
}
double y1(double y){
        return 0;
}
double y2(double y){
        return 1;
}
double rho(double x, double y){
        return 0;
}

int main(){

        int i,n=100;
        double xa=0, ya=0, h=1.0/n;

        double **u = new double *[n];
        for(i=0;i<n;i++) u[i] = new double [n];

        poisson2d( n, u, xa, ya, h, 0.001, x1, x2, y1, y2, rho);
        return 0;
}
\end{lstlisting}
这里设置左边界为电势为$0$的极板，右边界为电势为$1$的极板，上下边界电势均匀变化，所以会得到匀强电场。
得到数据以后，用gnuplot作图
\begin{lstlisting}
set term png
set output "poisson2d.png"
unset border
set xtics ("0.0" 0, "0.2" 20, "0.4" 40, "0.6" 60, "0.8" 80, "1.0" 99)
set ytics ("0.0" 0, "0.2" 20, "0.4" 40, "0.6" 60, "0.8" 80, "1.0" 99)
plot "poisson2d.txt" matrix with image
\end{lstlisting}
就会得到二维电势图\ref{fig:poisson2d}。
\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{poisson2d}
\caption{二维泊松方程：匀强电场。左侧极板电势为0，右侧极板电势为1。}
\label{fig:poisson2d}
\end{figure}

\section{波动方程}
偏微分方程这一章的算法思路主要是参考\cite{Hjorth-Jensen}，但其他部分，包括所有叙述、代码、例子，都是我自己做的。
\subsection{1维波动方程}
一维波动方程（比如绳波）的方程为
\begin{eqnarray}
\frac{ \partial^2 }{ \partial x^2 } u  = \frac{ \partial^2 }{ \partial t^2 }u,
\label{eqn:wave1d}
\end{eqnarray}
如果初始时1m长的绳子静止，绳子一端做简谐振动，另一端固定在墙上，那么初始条件和边界条件分别为
\begin{eqnarray}
u(x,t=0) &=& 0, \\
u(0,t) &=& sin(t), \\
u(1,t) &=& 0.
\end{eqnarray}
还需要一个初始条件，描述初始时刻绳子上各个质点的运动速度，
\begin{equation}
\frac{\partial}{\partial t} u |_{t=0} = 0.
\end{equation}
如果把微分方程用有限差分代替，即
\begin{eqnarray}
\frac{ u_{i+1,j} + u_{i-1,j} - 2 u_{i,j} }{(\Delta x)^2} + O((\Delta x)^2)
= \frac{ u_{i,j+1} + u_{i,j-1} - 2 u_{i,j} }{(\Delta t)^2} + O((\Delta t)^2),
\end{eqnarray}
如果取$\alpha = \frac{(\Delta t)^2}{(\Delta x)^2}$，显式的递推式为
\begin{equation}
u_{i,j+1} = (2-2\alpha) u_{i,j} + \alpha u_{i+1,j} + \alpha u_{i-1,j} - u_{i,j-1}.
\end{equation}
只要知道$u_{*,j-1}$与$u_{*,j}$，就可以递推得出$u_{*,j+1}$。

第一步前推时，我们只知道$u_{*,0}$，不知道$u_{*,-1}$，那么可以利用初始条件
\begin{equation}
\frac{\partial}{\partial t} u |_{t=0} = 0.
\end{equation}
翻译成差分，可以得到
\begin{equation}
u_{i,-1} = u_{i,1},
\end{equation}
带入递推式，得到
\begin{equation}
u_{i,1} = (1-\alpha) u_{i,0} + \alpha \frac{u_{i+1,0} + u_{i-1,0}}{2}.
\end{equation}
这样就可以完成前推了。
我还不清楚这个显式格式的稳定性，所以谨慎一点，将$\Delta t$取得小一点看看。
下面的代码中是一个函数，在给定初始条件，边界条件下，计算一维波动方程(\ref{eqn:wave1d})，并将各个时间节点的所有$u_{*,j}$作为一行输出到给定文件"file"中。
\begin{lstlisting}
void wave1d(    double xa, int nx, double dx,
                double (*boundary1)(double t), double (*boundary2)(double t),
                int nt, double dt, double (* u_initial)(double x),
                string file ){

        int i,j;
        double alpha = dt*dt/dx/dx;

        double *u_old = new double [nx]; // u_{i-1}
        double *u = new double [nx]; // u_i
        double *u_new = new double [nx]; // u_{i+1}

        u_old[0] = boundary1(0);
        u_old[nx-1] = boundary2(0);
        for(i=1;i<nx-1;i++) u_old[i] = u_initial( xa + i*dx);

        u[0] = boundary1(dt);
        u[nx-1] = boundary2(dt);
        for(i=1;i<nx-1;i++) u[i] = (1-alpha)*u_old[i] + alpha/2*( u_old[i+1] + u_old[i-1] );

        ofstream fp(file);
        if(!fp){
                cout<<"error: failed to open "<<file<<endl;
                exit(1);
        }
        for(i=0;i<nx;i++)fp<<u_old[i]<<"\t"; fp<<endl;
        for(i=0;i<nx;i++)fp<<u[i]<<"\t"; fp<<endl;

        for(j=2;j<=nt;j++){
                u_new[0] = boundary1( j*dt);
                u_new[nx-1] = boundary2( j*dt);
                for(i=1;i<nx-1;i++)
                        u_new[i] = (2-2*alpha)*u[i] + alpha * u[i+1] + alpha* u[i-1] - u_old[i];
                for(i=0;i<nx;i++)fp<<u_new[i]<<"\t"; fp<<endl;

                for(i=0;i<nx;i++){
                        u_old[i] = u[i];
                        u[i] = u_new[i];
                }
        }
        fp.close();
}
\end{lstlisting}
我们可以设置初始时绳子静止，左端点开始做简谐振动，右端点固定。
那么简谐波会从左往右传播，在右端点处反弹，反弹的波从右往左传播，与从左往右的波相干涉，形成驻波。
生成数据文件的代码如下：
\begin{lstlisting}

double u_initial(double x){
        return 0;
}

double boundary1(double t){
        return sin(10*t);
}

double boundary2(double t){
        return 0;
}

int main(){

        int i,j,nx=300;
        double xa=0, xb=3;
        double dx = (xb-xa)/nx;
        int nt=500;
        double ta=0, tb=5;
        double dt = (tb-ta)/nt;
        double **matrix;

        string file = "wave1d.txt";

        wave1d( 0, nx, dx, boundary1, boundary2,
                nt, dt, u_initial, file);

        // transpose the data matrix
        matrix = new double *[nt+1];
        for(i=0;i<nt+1;i++) matrix[i] = new double [nx];

        ifstream fp1( file );
        for(i=0;i<nt+1;i++)
                for(j=0;j<nx;j++)
                        fp1>> matrix[i][j];
        fp1.close();

        ofstream fp2( file );
        for(j=0;j<nx;j++){
                fp2<< xa + j*dx <<"\t";
                for(i=0;i<nt+1;i++)
                        fp2<< matrix[i][j]<<"\t";
                fp2<<endl;
        }
        fp2.close();

        return 0;
}
\end{lstlisting}
代码中最后一段是将数据文件的行和列做转置，方便用gnuplot画成不同的桢，然后生成动图gif文件。
gnuplot代码如下
\begin{lstlisting}
set term gif animate delay 1
set output "wave1d.gif"
set xrange [0:3]
set yrange [-2:2.5]
set xlabel "x(m)"
set ylabel "y(m)"
do for[t=2:502:1]{
        plot "wave1d.txt" u 1:t w l lw 3 t sprintf("string wave: v=1m/s, t=%f s", (t-2)*0.01)
}
set output
\end{lstlisting}
生成的动图在pdf文件中不方便展示，将在课堂上播放给学生看。

\subsection{2维波动方程}
2维波动方程如下
\begin{eqnarray}
\frac{ \partial^2 }{ \partial x^2 } u + \frac{ \partial^2 }{ \partial y^2 } u  = \frac{ \partial^2 }{ \partial t^2 }u,
\label{eqn:wave2d}
\end{eqnarray}
初始条件和边界条件为
\begin{eqnarray}
u(x,y,t=0) &=& g(x,y), \\
\frac{\partial}{\partial t} u(x,y,t) |_{t=0} &=& 0, \\
u(x=x_a, y, t) &=& Bxa(y,t), \\
u(x=x_b, y, t) &=& Bxb(y,t), \\
u(x, y=y_a, t) &=& Bya(x,t), \\
u(x, y=y_b, t) &=& Byb(x,t).
\end{eqnarray}
如果取$\Delta x = \Delta y = h$，则用有限差分代替2阶偏微分，可以得到
\begin{eqnarray}
&&\frac{ u^k_{i+1,j} + u^k_{i-1,j} - 2 u^k_{i,j} }{h^2} + O(h^2)
+ \frac{ u^k_{i,j+1} + u^k_{i,j-1} - 2 u^k_{i,j} }{h^2} + O(h^2) \nonumber\\
&=& \frac{ u^{k+1}_{i,j} + u^{k-1}_{i,j} - 2 u^k_{i,j} }{(\Delta t)^2} + O((\Delta t)^2),
\end{eqnarray}
其中$u^k_{i,j} = u(x_a + i*h, y_a + j*h, t_a + k*\Delta t)$。
记$\alpha = (\Delta t)^2 / h^2$，得到显式递推式
\begin{equation}
u^{k+1}_{i,j} = (2-4\alpha)u^k_{i,j} - u^{k-1}_{i,j} + \alpha ( u^k_{i-1,j} + u^k_{i+1,j} + u^k_{i,j-1} + u^k_{i,j+1} ).
\end{equation}
在初始时，利用初始条件$\frac{ \partial }{ \partial t} u |_{t=0} = 0$，可以得到 $u^{-1}_{i,j} = u^1_{i,j}$，得到递推公式的第一步
\begin{equation}
u^1_{i,j} = (1-2\alpha)u^0_{i,j} + \frac{\alpha}{2} ( u^0_{i-1,j} + u^0_{i+1,j} + u^0_{i,j-1} + u^0_{i,j+1} ).
\end{equation}
所以很容易将前面的1维波动代码改写为2维波动代码。
下面的代码实现二维平面波干涉，
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<fstream>
#include<cmath>
#include<string>

double source(double t){
        return sin(2*M_PI*t);
}

void wave2d(    double xa, double ya, int n, double h,
                double (*Bxa)( double y, double t ), double (*Bxb)( double y, double t ),
                double (*Bya)( double x, double t ), double (*Byb)( double x, double t ),
                int nt, double dt, double (* u_initial)(double x, double y),
                string file ){

        int i,j,k;
        double alpha = dt*dt/h/h;
        int ct = n/2;

        ofstream gnu("wave2d_interference.gnu");
        gnu<<"set term gif animate delay 5\n";
        gnu<<"set output \"wave2d_interference.gif\"\n";
        gnu<<"set xrange [0:10]\n";
        gnu<<"set yrange [0:10]\n";
        gnu<<"set zrange [-0.5:1]\n";
        gnu<<"unset border\n";
        gnu<<"unset xtics\n";
        gnu<<"unset ytics\n";
        gnu<<"unset ztics\n";
        gnu<<"set view 45, 30, 1, 1\n";

        double **u_old = new double *[n];
        for(i=0;i<n;i++) u_old[i] = new double [n];
        double **u = new double *[n];
        for(i=0;i<n;i++) u[i] = new double [n];
        double **u_new = new double *[n];
        for(i=0;i<n;i++) u_new[i] = new double [n];

        // calculate u^0_{i,j}
        for(i=0;i<n;i++){
                u_old[0][i] = Bxa( ya + i*h , 0 );
                u_old[n-1][i] = Bxb( ya + i*h, 0 );
                u_old[i][0] = Bya( xa + i*h, 0 );
                u_old[i][n-1] = Byb( xa + i*h, 0 );
        }
        u_old[ct-15][ct] = source(0);
        u_old[ct+15][ct] = source(0);
        for(i=1;i<n-1;i++){
                for(j=1;j<n-1;j++){
                        if((i!=ct-15 || j!=ct) && (i!=ct+15 || j!=ct) )
                        u_old[i][j] = u_initial( xa + i*h, ya + i*h );
                }
        }
        // output u^0_{i,j}
        string file_frame = file + "000.txt";
        ofstream fp(file_frame);
        if(!fp){
                cout<<"error: failed to open "<<file_frame<<endl;
                exit(1);
        }
        for(i=0;i<n;i++){
                for(j=0;j<n;j++){
                        fp<< xa + i*h <<"\t"<< ya + j*h<<"\t"<< u_old[i][j]<<endl;
                }
                fp<<endl;
        }
        fp.close();

        gnu<<"splot \""<<file_frame<<"\" with pm3d t sprintf(\"t=%f s\", 0)\n";


        // calculate u^1_{i,j}
        for(i=0;i<n;i++){
                u[0][i] = Bxa( ya + i*h , dt );
                u[n-1][i] = Bxb( ya + i*h, dt );
                u[i][0] = Bya( xa + i*h, dt );
                u[i][n-1] = Byb( xa + i*h, dt );
        }
        u[ct-15][ct] = source(dt);
        u[ct+15][ct] = source(dt);
        for(i=1;i<n-1;i++){
                for(j=1;j<n-1;j++){
                        if((i!=ct-15 || j!=ct) && (i!=ct+15 || j!=ct) )
                        u[i][j] = (1-2*alpha)*u_old[i][j] + alpha/2*( u_old[i-1][j] + u_old[i+1][j] + u_old[i][j+1] + u_old[i][j-1] );
                }
        }

        // output u^1_{i,j}
        file_frame = file + "001.txt";
        fp.open(file_frame);
        if(!fp){
                cout<<"error: failed to open "<<file_frame<<endl;
                exit(1);
        }
        for(i=0;i<n;i++){
                for(j=0;j<n;j++){
                        fp<< xa + i*h <<"\t"<< ya + j*h<<"\t"<< u[i][j]<<endl;
                }
                fp<<endl;
        }
        fp.close();

        gnu<<"splot \""<<file_frame<<"\" with pm3d t sprintf(\"t=%f s\", "<<dt<<")\n";

        for(k=2;k<=nt;k++){
                // calculate u^{k}_{i,j}
                for(i=0;i<n;i++){
                        u_new[0][i] = Bxa( ya + i*h , k*dt );
                        u_new[n-1][i] = Bxb( ya + i*h, k*dt );
                        u_new[i][0] = Bya( xa + i*h, k*dt );
                        u_new[i][n-1] = Byb( xa + i*h, k*dt );
                }
                u_new[ct-15][ct] = source(k*dt);
                u_new[ct+15][ct] = source(k*dt);
                for(i=1;i<n-1;i++){
                        for(j=1;j<n-1;j++){
                                if((i!=ct-15 || j!=ct) && (i!=ct+15 || j!=ct) )
                                u_new[i][j] = (2-4*alpha)*u[i][j] - u_old[i][j] + alpha*( u[i-1][j] + u[i+1][j] + u[i][j+1] + u[i][j-1] );
                        }
                }
                for(i=0;i<n;i++){
                        for(j=0;j<n;j++){
                                u_old[i][j] = u[i][j];
                                u[i][j] = u_new[i][j];
                        }
                }
                // output u^k_{i,j}
                file_frame = file + (char)(48+(k/100)%10) + (char)(48+ (k/10)%10) + (char)(48+ (k%10)) + ".txt";
                fp.open(file_frame);
                if(!fp){
                        cout<<"error: failed to open "<<file_frame<<endl;
                        exit(1);
                }
                for(i=0;i<n;i++){
                        for(j=0;j<n;j++){
                                fp<< xa + i*h <<"\t"<< ya + j*h<<"\t"<< u_new[i][j]<<endl;
                        }
                        fp<<endl;
                }
                fp.close();
                gnu<<"splot \""<<file_frame<<"\" with pm3d t sprintf(\"t=%f s\", "<< k*dt <<")\n";
        }
        gnu<<"set output"<<endl;
        gnu.close();
}

double u_initial(double x, double y){
        return 0;
}

double Bxa(double y,double t){
        return 0;
}
double Bxb(double y, double t){
        return 0;
}
double Bya(double x, double t){
        return 0;
}
double Byb(double x, double t){
        return 0;
}

int main(){

        int i,j,n=300;
        double xa=0, xb=10;
        double ya=xa, yb=xb;
        double h = (xb-xa)/n;
        int nt=300;
        double ta=0, tb=5;
        double dt = (tb-ta)/nt;
        cout<<"dt="<<dt<<endl;

        string file = "wave2d/interference";
        wave2d( xa, ya, n, h,
                Bxa, Bxb, Bya, Byb,
                nt, dt, u_initial, file );

        return 0;
}
\end{lstlisting}
这个代码既做微分方程的计算，生成数据文件，也生成gnuplot脚本，执行完这个程序以后，再运行生成的gnuplot脚本，就会生成2维平面波干涉的gif动图。
